Preprint, March 1995 


Threshold dynamics and strength-duration curves 


Roberto Sudrez-Antola 
Facultad de Ingenieria, Universidad Catélica del Uruguay and Direccién Nacional de 
Tecnologia Nuclear, Ministerio de Industria, Energia y Mineria, Montevideo, Uruguay 


Abstract: The mathematical modelling of electric current flow in excitable cells (fibres 
and interrelated cells networks) is briefly reviewed and three fundamental problems of 
the stimulation by external electrodes are defined and posed, in terms of suitably space- 
averaged fields. Then, it is given an outline of a new procedure to study just-threshold 
dynamics when certain target elements in excitable tissues are stimulated by external 
electrodes. The concept of “region of influence” of an external field of applied electric 
current over the target element is introduced as a first step towards a nonlinear modal 
analysis of membrane’s polarization. In a suitable defined mode amplitude space, a 
threshold manifold, a decaying set and an amplifying set are introduced to study the 
emergence of an action potential. It is then outlined a general method to derive strength- 
duration formulae, giving the derivation of Lapicque-Hill’s formula as an example. 


A. Introduction. 


Strength duration curves measure the excitability of a certain target element in a given 
electrode-tissue system. The target element may be a single cell (such as a nerve or 
skeletal muscle fibre) or even a whole network of electrically coupled cells (such as the 
myocardium or visceral smooth muscle). The response of the target element is the all or 
none emergence of a propagated action potential. The results of the experiments show 
that the parameters of the strength-duration curve depend not only of the excitability 
properties of the membranes of the target element, but depend also of certain geometric 
and bioelectric parameters of the system formed by the electrodes, the tissues and the 
target element. For example: radius of the excitable fibre, distance between the electrode 
and the target element, size and shape of the active electrode (if it is located close 
enough to the target element), amongst others (1, 2, 3). So, strength duration curves 
measure something that may be called a “global excitability” of the target element in the 
framework of a given electrode-tissue system, by opposition to the “local excitability” of 
the membranes of the target element. This local excitability could be measured, at least 
in principle, determining the strength-duration curves for the uniformly polarized 
membranes (by a patch-clamp technique, for example) and relating its parameters with 
the known properties of the ionic channels. The explanation of the experimental results 
that show a relation between the phenomenological parameters of the strength-duration 
curve and certain parameters of the electrode-tissue system has remained as an open 
problem from more than seventy years. This problem has more than an academic 
interest, because it is closely related with two basic questions of clinical electro 
stimulation: Which elements will be excited by a given current pulse in a given 
electrode-tissue system? How to reach a selective stimulation of certain target elements 
in a whole excitable tissue? In order to get some kind of solution to the aforementioned 
problem, it seems necessary to solve first the following three problems of the electric 
stimulation by external electrodes: 


(1) Given a certain electrode array in a certain volume conductor formed by 
biological tissues, to find the resulting field of electric current density. 


Preprint, March 1995 


(2) To obtain a relation between this current density field and the space-time 
pattern of polarization of excitable membranes, including the field of 
transmembrane voltage Vy and the whole set of activation, recovery and 
adaptation variables of the ionic channels, represented by the vector field w. 


(3) To identify certain spatial patterns of (Vy,,w) that in some sense could be 
considered as “threshold states” of the excitable membranes in the target 
elements. 


B. The electric current flow in excitable tissues. 


Let us first consider the electric current flow produced by the external electrodes. If the 
emergence of an action potential were a short range spatial phenomenon, we would need 
a fairly detailed description of the flow in the environment and in the interior of the 
excitable cells. But, fortunately, this emergence is a long range spatial phenomenon 
under common physiological conditions. So, it is possible to work with suitable space- 
averaged fields, smoothing unnecessary details and greatly simplifying the mathematical 
formulation of the three problems of electro stimulation by the use of an extension of 
the so-called bidomain model. 

The bidomain model is a phenomenological approximation (4) originally devised to 
describe the electric current flow in cardiac tissue in the framework of 
electrocardiography and magnetocardiography, but it can be restated to study threshold 
dynamics in any electrode-tissue system (5, 6). 

It assumes that the discrete nature of the tissues can be ignored in favour of a simplified 
continuum approach. 

Three continua are postulated in the original version of the theory, each one filling the 
same spatial region: an intracellular domain an interstitial domain and the continuum of 
excitable membranes. In each continuum, several fields are introduced that verify a 
coupled set of nonlinear partial differential equations. 

To study the problems of excitability it is convenient to substitute the intracellular and 
the interstitial domains by a single ohmic volume conductor, the equivalent volume 
conductor. As this new continuum is independent from the continuum of excitable 
membranes, it is possible to pose and solve the problem of finding the electric current 
distribution produced by an array of external electrodes as in any other ohmic volume 
conductor. 

Now we must turn to consider the target element. If it is a single excitable fibre, it can 
be represented by a global, spatially one-dimensional and nonlinear cable model (7). 
The cable may be idealized by a certain curve in the region occupied by the volume 
conductor. If s represents the arc length along that curve, a fairly general version of the 
nonlinear cable equations for a uniform fibre is the following (6): 


OV = oN 0°O@ 
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Here + aw = Bu ou is membrane’s linear time constant, Ry is the resistance of the unit 


area of membrane, Cy is the capacitance per unit area, J Vy. W) is the ionic current 
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density through membrane’s ionic channels, Aq is membrane’s space constant, 
F (V,,,w) represents the kinetics of the activation and recovery variables of the ionic 
channels and its scalar components are given by expressions of the Hodgkin and 
Huxley’s type (but that vary from one membrane to another), and t represents time. 

The coupling between the state variables (V,,,w) of the nonlinear cable equations 


[la], [1b] and the applied electric current field in the equivalent volume conductor is 
2 


2 
given by the source term? y” — where Z is equivalent to the activating 


sg” 6s" 
function introduced by Rattay (8). 
2 


Under quasi-steady conditions ge As F(s) I(t), where I(t) is the total applied current 


and F(s) is what we call the form factor (6). When the fibre is bounded, another 
coupling with the external current field appears through the boundary conditions at the 
end points of the fibre (6). 


C. Modal analysis of threshold dynamics: a brief outline. 


If we want to continue with an analytical approach to threshold dynamics, we need a 
new idea that allows us to drastically simplify the original nonlinear mathematical 
problem. This idea happens to be the region of influence of the applied current field 
over the membranes of the target element. 

If the target element is a single fibre, and if there is a single active electrode (cathode or 
anode), close enough to the fibre, this region is simply the interval of influence of the 
electrode over the fibre (6). 

If the target element is a functional syncytium and if there is a single active electrode, 
like in cardiac pacing, we have a whole three-dimensional region of influence of the 
electrode over the tissue (6). 

Using the region of influence it is possible to pose and then solve a spatial Sturm- 
Liouville problem (6) (with Dirichlet’s boundary conditions for a fibre and with mixed 
Dirichlet’s and Neumann’s boundary conditions for a bounded continuum of excitable 
membranes). 

The fields Vy and wcan be developed as a series of functions in terms of the 
corresponding eigenfunctions or modes, with mode amplitudes that depend only of time. 
If we approximate membrane’s ionic current density J(V,,,w) and the scalar 
components of the function F (V,,,w) by suitable polynomials in terms of the deviations 
of the state variables from the rest state, the nonlinear partial differential equations [1] 
can be substituted by an infinite set of ordinary nonlinear differential equations in terms 
of the amplitudes of the modes. 

Moreover, it is possible to truncate these equations for the mode amplitudes from a 
certain mode index onwards. 

In that case the whole spatial distribution of the state variables Vm and w of the target 
element membranes can be represented by a single point in a suitable finite-dimensional 
space of mode amplitudes. The rest state is represented by the origin. 

Furthermore, the widely separated time scales of evolution of activation variables, of 
membrane potential and of recovery variables, allows the partition of that space in two 
sets: a decaying set and an amplifying set. 

The rest state is always in the decaying set. 


Preprint, March 1995 


If at the end of an applied pulse the representative point is in the decaying set, the 
excitation will fail and the membranes of the target element will return to rest. 

But if it is in the amplifying set, an action potential will be produced. 

The amplifying set is separated from the decaying set by a more or less well defined 
transition region that can be idealized by a certain threshold manifold in the space of 
mode amplitudes. 

Each point of this manifold corresponds to a threshold spatial distribution of the 
variables Vy andw, so that we obtain a definite answer to the third posed problem of 
the excitation by external electrodes. 

Moreover, to each history of applied current corresponds a certain orbit in the space of 
mode amplitudes (the effect of the external current on the rate equation for a given mode 
amplitude appears as a source term proportional to I (t)). 

This allows us to solve the second problem of electro stimulation. Considering the 
intersection between the decaying set and the region reachable from the origin by 
cathodic or anodic pulses, applied by a single active electrode, phenomena like cathodic 
make excitation, anodic make excitation, anodic block from the cathode, anodic break 
excitation and cathodic break excitation can be analysed by a thorough examination of 
the location of the reachable subset of the threshold manifold in mode amplitude space. 
A rectangular pulse of externally applied current will be just-threshold if and only if, at 
the end of the pulse the point which represents the state of the membranes of the target 
element belongs to the threshold manifold. 

Now, linearizing the dynamics of the mode amplitudes in the decaying set, simplifying 
(if necessary) the threshold manifold and taking into account the aforementioned just- 
threshold condition for a rectangular pulse of applied electric current, it should be fairly 
straight-forward to derive analytical expressions for the strength-duration curves, taking 
into account the effect of any number of activation and recovery variables as needed. 


D. A derivation of Lapicque-Hill’s formula for cathodic make excitation. 


The simplest example of the general procedure outlined in section C is afforded by a 
straight and uniform fiber whose membrane can be considered with its recovery 
variables frozen at its rest values and with its activation variables always relaxed to its 
steady-state values. 
Then, the ionic current density of the membrane is a function of Vy alone and can be 
fairly well approximated by a cubic polynomial. In that case we can work with a single 
nonlinear partial differential equation (Nagumo-Rattay’s equation) (6). 
If we introduce the further restriction of considering only cathodic make excitation by a 
small convex electrode, then the amplitude A, of the first mode happens to be dominant, 
and as first approximation it can be uncoupled from the other mode amplitudes (6). 
With this simplification we obtain a one-dimensional mode amplitude space (the line 
that corresponds to A;) in which the dynamics is given by the following equation (6): 
2 
Gy 5 Gee a Ana 
dt id 
Here ? is the length of the interval of influence of the electrode over the fibre. 
It is nearly equal to 6z (6), where z is the distance between the centre of the electrode 
and the axis of the fibre. 
F, is the projection of the form factor F(x) onto the first eigenfunction. It is inversely 


proportional to @° (6). 


yA, +b 9, A, -69,, A) tay EI [2] 
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The positive parameters b and c stem from the cubic polynomial that gives the ionic 
current density through fibber’s membrane, while 911; = 1,2 and 9411; = 1,5. 

If we put /(t) equal to zero, the resulting equation has three equilibrium points: the 
origin (stable), a threshold point Aj, (unstable) between the other two, and an excited 
state Aj. (stable). In this case the decaying set is given by A; < Aju, the amplifying set is 
given by A; > Aj, and the threshold manifold is given by 


2 
A = PH (PS) 1 fy, aw 3] 
7 2cH 1, 2¢4K 11, CHA, c 
Linearizing the equation [2] in the decaying set we obtain: 
dA, n° Ae 2 
ies 1+ zB A, +Ay F I(t) [4] 


Now, consider a rectangular pulse of amplitude / and duration d that begins with the 
membrane at rest (that is A;(O) = 0). 


The condition for a just-threshold pulse is that: Ai (d) = Atu [5] 
Integrating [4] for a rectangular current pulse and taking into account the just-threshold 
. R 
condition [5] we finally obtain Lapicque-Hill’s formula: eer [6] 
f - | 
2 2 
A 

The rheobase is: R=|14+7 a 7: [7] 

CO Ay? F 

are e 
The system’s time constant is: i= — [8] 
TE Ans 
1+ ry 


For a convex electrode: ¢ = 6(r) +h) being ro the electrode’s radius and h the distance 


between the electrode’s surface and the target fibre (6). 

It is possible to make a somewhat similar derivation of Lapicque-Hill’s formula for the 
case of cathodic make stimulation of myocardium by a convex electrode, but the 
mathematical analysis is less simple (6). 

The fibre model considered in this section can be improved to take into account that 
most fibres are wrapped by a purely resistive sheath of connective tissue. The 
importance of this sheath in the distribution of external current through fibre’s 
membrane grows as the electrode gets closer to the fibre. It is possible to derive 
Lapicque-Hill’s formula also in this case (6). 


E. Final comments and conclusions. 
Formulae [7] for the rheobase and [8] for the time constant solve the problem posed in 
the introduction, for a single uniform fibre excited by a convex cathode, when the 
activation variables can be considered relaxed to equilibrium with membrane potential 


and the recovery variables can be considered frozen. 


5 
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Even when the improved formulae, which take into account the effect of a resistive 
sheath, seem to agree well with the available experimental data, a proper test of the 
predictions of the theory can be made only through thoroughly designed new 
experiments. 

Some preliminary work has been done by Michel Borde and Sebastian Curti, from the 
Department of Neurophysiology, Faculty of Medicine, University of the Republic, 
Uruguay. 

The results show a fairly good agreement with the theoretical predictions when the 
rheobase, time constant and limit threshold charge are measured in vitro against the 
distance between a point electrode and the median giant nerve fibre of the earthworm 
(when completed, they will be reported elsewhere). 

The method outlined here can be applied, in principle, to the excitation of nerve, skeletal 
muscle, myocardium and smooth muscle, not only when the active electrode is convex 
but also when it is concave. 

The main shortcoming of the proposed general procedure is that it is unable to cope 
with the growth of the action potential, being restricted to just-threshold dynamics. 
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FROM SYNERGETICS 
TO BIOLOGICAL EXCITABILITY 


Anibal C. Sicardi-Schifino and Roberto E. Suadrez-Antola 


1 - Introduction 


The emergence of a propagated disturbance (action potential) in a long nerve fiber 
or in the myocardium, when this tissue is electrically stimulated by an external 
electrode, is an essentially nonlinear and localized (space-dependent) phenomenon. 
On the contrary, when the tissue is stimulated by microelectrodes in laboratory 
experiments, one can produce an uniform polarization. In this case the phenome- 
non is nonlinear but nonlocalized and it exists a theory based in the Fitzhugh- 
Nagumo equations, which are a simplified version of Hodgkin and Huxley ones 
(Jack et al. 1983; Volkenstein 1983). These equations are of reaction-diffusion type 
and present well known nonlinear effects. In this case, the propagation of the 
action potential can be understood in terms of solitary-like waves. When the 
stimulation is time periodic, Hopf bifurcations and chaotic behavior appear. 
However, the electrical stimulation of excitable tissues by external electrodes is 
fairly important from a medical and electrophysiological standpoint. (For example, 
functional electric stimulation of nerve and skeletal muscle fibers and chronic 
electric pacing of the heart and parts of the nervous system.) In this kind of 
stimulation and in the synaptic one, membrane polarization is essentially localized. 
The classical concept of a uniform threshold potential can’t be applied to study the 
emergence of a propagated disturbance (action potential) in the tissue. Of course, 
as in the uniform case (space clamped), this emergence is essentially nonlinear. 
But now, the success or failure of the excitation of a given target element depends 
not only of the history of global applied current, but of the position of that element 
in the electric current field produced by the electrodes. As a consequence of this 
spatiotemporal dependence, two closely related problems arise: to establish which 
elements are excited in a given electrode-tissue system and how to reach a selective 
stimulation of certain target elements using the spatiotemporal properties of 
electrode-tissue systems. Several numerical simulations have been done to study 
the emergence of an action potential for nerve fibers and other tissues stimulated 
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by external electrodes. Recently, using the tools of Synergetics, we began to develop 
{in a series of papers: Suarez & Sicardi 1996; Suarez et al. 1997) a nonlinear 
analytical theory for this spatially localized phenomenon, which is reviewed in this 
lecture (while the biological results obtained are reviewed in other paper [Suarez & 
Sicardi 1997] in this same volume). The goal of our work is to contribute to the 
establishment of a bridge between the local biophysical approach to excitability 
and the global approach of clinical electrophysiology. Until now, we only have 
studied the threshold conditions for the generation of action potential and not the 
propagation of a finite disturbance (that we will consider in a future work). 
Moreover, here we only consider the simple case of a single unmyelinated nerve (or 
muscle) fiber and assume that fast membrane activation variables are relaxed to 
equilibrium. Previously to present our theory, we review Haken's slaving principle 
and abridge the main ideas of Synergetics. (Haken 1978; Haken 1983) 


2 - Mathematical tools: synergetics 


Synergetics is the study of self-organization phenomena and pattern formation. 
That is, it studies the emergence of dissipative structures in nonequilibrium open 
systems. In this phenomena, the cooperation of many individual (microscopic) 
parts of a system leads to dramatic changes on the global (macroscopic) structure 
of the system, in a self-organized manner. The typical example is the Bénard’s 
convection instability in fluids: We consider a fluid layer heated from below (and 
kept at a fixed temperature from above). For small temperature gradients, heat is 
transported by conduction (from the lower to the upper surface) and no macrosco- 
pic fluid’s motion shows up. When the temperature gradient exceeds a certain 
critical value, suddenly a macroscopic fluid’s motion (the convection) occurs, in the 
form of rolls. This phenomenon is an example of bifurcation, where the temperature 
gradient is the control parameter. (Haken 1978; Guckenheimer & Holmes 1983) 

In principle, in order to describe a pattern formation, as in the case of convec- 
tion instability, we need to solve a system of nonlinear partial differential equations 
{the Navier-Stokes equations, in the considered example) with boundary conditions. 
But, making a nonlinear modal analysis (i.e. a Fourier transform), we can convert 
them into infinite coupled (nonlinear) ordinary differential equations. For this, first 
we must state a Sturm-Liouville problem, that is to obtain an appropriate 
orthonormal set of eigenfunctions y,(x) which satisfy the boundary conditions. 
Then, we must expand all the physical variables of the problem in terms of the 
eigenfunctions y (x) (ie. in the form ¥_4(t)w,(x) where the coefficients of the 
expansion A, are the mode amplitudes.) Thus. projecting the partial differential 
equations onto the eigenfunctions (that is, performing the corresponding functional 
scalar product). we get a system of infinite nonlinear ordinary differential equations 
for the mode amplitudes 4 ft) 
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Generally. only a few modes. called order parameters. are relevant. They 
dominate the system dynamics. While the other modes are slaved by them. Thus. 
using the Haken’s slaving principle, it is possible to eliminate the slaved modes, 
reducing the problem to solve a few (coupled) nonlinear ordinary differential 
equations (the Lorenz equations. for the case of convection instability). These 
resulting equations can show a very complex behavior. with bifurcations and 
(sometimes) chaos. For example, as it is well known, the Lorenz equations have 
several bifurcations and, for certain values of control parameters, their solutions 
are chaotic. In order to show how the slaving principle works, we consider the 
following two-modes example: 


Dia A, 4A, , Beard, tba, 


where t, is small and its sign is indefinite, while t, is always positive and big, so 
that, for long times, it is approximately d4,/dt =0. Therefore, we may write 


and 
6 
=-1,4, —« (2.1) 


In this form 4, is slaved by A, and it is eliminated. Now, we must solve only the 
equation for the dominant mode (or order parameter) 4,, i.e. Eq. (2.1). If we assume 
4/,=1 and t,=-a, we obtain the normal form of the pitchfork bifurcation (Haken 
1978; Drazin 1992), 

Fnas- 4. (2.2) 

For & negative, the only (stable) equilibrium point is 4,=0 while for & positive 
there are three equilibrium points: one unstable (4,=0) and two new (stable) 
solutions, 4,=¥a and 4,=-¥Ja. In the critical value a=0, the pitchfork bifurcation 
occurs (Q is the control parameter). An example of pitchfork bifurcation is the 
Benard instability, mentioned before. In this case the equilibrium 4, =0 represents 
the rest state of the fluid (whit only heat conduction), while the other two equilibria 
(4,=Va and 4,=-Va ) represent the convection rolls. 

In the above example of bifurcation (the pitchfork bifurcation), the control 
Parameter was multiplicative. Now we will consider an example with an additive 
control parameter and so, we will generalize the slaving principle for this case. 

We begin considering a system described by the following equations: 


dg q, 44, 
HUAI «GEFs» =H Fah « as 
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where F is the the control parameter, which is always positive. The trivial 

equilibrium is g,=9,=0, g,=F. which is unstable if F > | (but it is stable for F 
between 0 and 1). Making the change of variables 

q,*+9, q,-49, 

Las. 4S 

_ =. 

we obtain the nonlinear but homogeneous system 


» 4=9,-F. 


aA dA dA A 
vals 44, 5 eile te a —— = =. (2.4) 
In this form, we transform the additive control parameter F in a multiplicative 
one, (F—1). (We note that, for long times, it is approximately dd, /dt = d4,/dt = 0). 
Now, using the slaving principle, we may eliminate the variables 4, and 4,, which 
are slaved by 4,. Therefore, it is approximately 4,=0 and 4, =-.4)/2. Then, we get 
the following equation for 4,, which has the form of Eq. (2.2) and so represents a 
pitchfork bifurcation, 


Sha (F-1)4,-440. 25) 


The same method may be used in other examples with an additive control 
parameter (that is, in cases of inhomogeneous nonlinear equations, which may be 
transformed in homogeneous ones, making a change of variables that eliminates 
the inhomogeneous terms). In all the cases, we reduce the problem to one with a 
multiplicative control parameter. 


3 - General equations 


We begin modeling a usual experimental situation, the stimulation by external 
electrodes of a long nerve fiber. In a first approach to the problem, we can neglect 
the activation variables effects and myelin’s ones. Thus, to sketch the phenomenon. 
we consider a long unmyelinated nerve fiber with a stimulating active electrode (a 
cathode or an anode) placed relatively near it. And we assume that membrane’s 
fast activation variables are relaxed to their equilibrium values with membrane 
potential. Using the charge conservation and Ohm's law for the internal potential 
¥ =V_,+(x,t)+¢, and assuming the most simple approximation for the membrane 
resistance (considering only one recovery variable), we obtain the following 
equations, (Suarez & Sicardi 1996; Jack et al 1983; Volkenstein 1983) that 
describe the polarizing effect of the external electrode (depolarizing effect if the 
electrode is a cathode), 


é 
=k —+f(v)-aws+i_— (3.1) 
ex 
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cw 
— =Biv-y) (3.2) 
cot 


These are the usual Fitzhugh-Nagumo equations, (Britton 1986: Tuckwell 1988) 
but with the additional (Rattay's [Rattay 1987]) term 22 (°4/ax*) that take account 
localized polarization effects of the external potential $(x,r). We can call Fitzhugh- 
Nagumo-Rattay equations to Eqs. (3.1)-(3.2). Here, v(x,r) is the difference between 
membrane voltage V_ and its resting value V_,, while 4, (fiber’s linear space 
constant) and t_ (membrane linear time constant) are characteristic length and 
time, respectively. Moreover, the recovery variable w(x,1) is also referred to its rest 
value, while we assume that rest of the membrane’s ionic current density is 
approximated by a cubic polynomial, so that 


fv) =-v4+bv' -ev' (3.3) 


where 6, c, &@, B and y are positive constants. 

Assuming that the external volume conductor in which the target fiber is 
immersed is purely resistive, with a conductivity G, a quasistatic approximation 
works. Thus, it is always possible to factorize the Rattay’s activating function 
0°9/&? as 


a 
—=F(x)l(t) (3.4) 
& 


where /(t) is the total injected external current and F(x) is a space dependent 
function, that we call in the following, the form factor. We note that this form factor 
depends on the space distribution of the electric current density and the target 
fiber’s position in this field. In the case of a point cathode, for example, we obtain 


4 (2x'-z') 
F(x)=(4nG) ——> (3.5) 
(x +z) 

where z is the distance between the point electrode and the fiber, while x is the 
Position along the fiber, being x=0 the fiber point nearest to the cathode. Ob- 
viously, in this case, /(¢) is always negative because it is a cathodic current. Then. 
the sign of F(x) and the activating function's one are opposites. 

In our first paper (Suarez & Sicardi 1996), which we review in Section 4, we 
have simplified the Fitzhugh-Nagumo-Rattay equations (3.1)-(3.2), assuming that 
the slow recovery variables are frozen in their resting values. That is, we assume. in 
this first approach, that w(x,r)=0. Therefore, in this case, the system (3.1)-(3.2) 
reduces to 


x 


2 ay 2 
tT a ah as (3.6) 
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However, in a further approach, (Suarez et al 1997) which we review in Section 
5, we consider the full system (3.1)-(3.2). Now, in both cases, we must make a 
nonlinear modal analysis of the threshold dynamics, given by Eq. (3.6) (in the 
simplified approach of Section 4) or by Eqs. (3.1)-(3.2) (in the full approach of 
Section 5). 

But, in order to carry out a nonlinear modal analysis of the threshold dynamics, 
we need to construct an appropriate orthonormal set of eigenfunctions. That is, we 
must state a Sturm-Liouville problem 


ie v=uv (3.7) 


with boundary conditions at the end points of a finite space interval. However, it is 
inadequate to identify this interval with the fiber’s length, because the fiber’s 
depolarization zone (due to the electrode’s action) in the threshold is usually fairly 
smaller than the fiber itself. Then, first at all. we must find the appropriate fiber’s 
space interval to state the boundary conditions. Obviously, this interval should be 
directly related with the localized depolarization of the target fiber. So, the most 
correct space interval to use would be the electrode’s influence interval, that is the 
length of the fiber's depolarization zone, at whose end points the membrane 
potential should be zero. But this length isn't well defined. However, as we will see 
below, it is always possible to estimate the length / of this influence interval. Thus, 
we can use this estimated influence interval, putting the membrane voltage equal 
to zero in its end points. Therefore, we can formulate a Dirichlet problem for 
equation (3.7) and, taking the coordinate origin at the middle of the influence 
interval, we get the following eigenvalues and eigenfunctions 
2.2.42 
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The boundary conditions that we state correspond to threshold cases (that is. to 
study the emergence of the action potential). These conditions must be modified in 
order to study the propagation of a finite action potential. We will consider these 
modifications in a future work. 

We note that the first mode (n=1) corresponds. in the cathodic case. to a pure 
depolarizing mode in the influence interval. While in the anodic case. it corres- 
ponds to a pure hyperpolarizing mode in this interval. That is, they are typical 
distributions of membrane polarization along the fiber axis. while the other modes 
(n=2.3....) represent the detailed corrections to these profiles. essentially the lateral 


v,(2)= sin when » is even (3.8) 
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“wings” of opposite polarization. For these reasons, we can hope that the solution of 
the complete Fitzhugh-Nagumo equations (3.1-3.2) be expressible as a linear 
combination of these biophysically significant modes. Moreover, we can expect that 
mode one’s amplitude be significantly greater than the other ones. 

Then, if we introduce the ansatz 


(x)= DA (oy, (x), 
w(x.t)= DB (t)y,(x) . 


in Eqs. (3.1-3.2), we obtain, for each mode n 
2 2,2 


dA z 
a ll 7 *)4,+609 44 -cD9 AAA —aB, +2 1(t)F. (3.9) 
?. 


weqtr 
a ear 
dB 
©. 4,-10, (3.10) 
where 
3, = [iv,@v, (Od. 9. = fy Wo, cv, (Oy, (de 
and 


= Fy, CO F ede . (3.11) 


In each mode-equation (3.9), MIE, is a source term proportional to applied 
current and to the form factor’s projection onto corresponding mode eigenfunction. 
We note that these source terms act directly on the membrane potential modes, but 
only indirectly on recovery variable modes. Moreover, we note that the product Br, 
is a small dimensionless number. 

Now, we consider how to estimate the length / of the influence interval in the 
cathodic case. The anodic case is analogous. The cathode’s influence interval is 
determined by the form factor F(x), defined in Eq. (3.4). In the point cathode’s case 
this form factor is given in Eq. (3.5) and its dependence with x is shown in Figure 1. 
For any finite symmetric convex electrode, the behavior of the form factor F(x) is 
qualitatively like point electrode’s one. That is, it has the same general traits shown 
in Figure 1. Moreover, when |x| is enough large, the form factor F(x) of any finite 
electrode is approximately equal to the point electrode’s one, given in Eq. (3.5) 
Thus, let us begin considering a point electrode (see Figure 1). In the interval 
~*, SxSx, the form factor F(x) is negative and so, leads to depolarize fiber’s 
membrane. For x < -x, and for x > x, there are two lateral “wings” where F(x) is 
Positive and then, tends to hyperpolarize the membrane. While the depolarization 
effect is well localized. the hyperpolarization effect in the “wings” is somewhat 
diffuse. As Figure 1 shows. each lateral “wing” presents a maximum [in +r,,) and 
an inflection point (in +x,). From x=+x, outward F(x) tends to zero as |x|". Then. a 
good estimation of the influence interval (see Figure 1) is the interval (-x.x) 
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= og, oa oe (3.12) 


dx 


All of this may be generalized to any bounded current sources distribution. 
However, let us consider always in the following symmetric form factors. so that 
F, =0 for any even n. 


Figure 1. The graphics of the Form Factor for a point electrede and the construction of the 
influence interval. 


4 - Threshold dynamics without recovery variables 


In our first paper (Suarez & Sicardi 1996), we freeze the recovery variable at its 
rest value, that is B =0 for every n. Now, we abridge the results of this paper. In 
this case, the system (3.1)-(3.2) reduces to Eq. (3.6). Then, Eqs. (3.9)-(3.10) reduce 
to Eq. (3.9) with 8. =0. As 4, is equal to zero for ¢=0, for all n, and the relaxation 
constant is minimum for n=1, if the only non negligible form factor is F,, then only 
the mode n=l is directly excited, meanwhile the other modes are slaved by this 
dominant mode, being 


ae 
_ - y" [9.4 +N EF] (4.1) 


A, =(I+ 
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for the modes with n > | (n=3,5,7,...).Then, Eq. (3.9) reduces in this case (at third 
order) to the following equation (for the mode one amplitude A} : 


£e hao 4) EMF, (4.2) 


where 
a a 
2 *)A,-b'A +c'A, (4.3) 


A 4,)=(1+ 


22,2 


ank 
b'=63,,.  c'=c9,,,,-26° D9? TT 


lea 


being the summation over the odds values of n, from 3 onward. Moreover, as was 
seen before, for biophysical reasons the slaved modes amplitudes for n#1 are small 
relatively to 4. This assumption can be further sustained by numerical analysis, 

5.00 


4.00 


3.00 


2.00 


0.00 -2.00 4.00 6.00 8.00 10.00 
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Figure 2. Bifurcation diagram for Eq. (4.2), showing a saddle-node bifurcation at the 
critical influence length. 
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for typical values of membrane parameters (for example. with the values given by 
Jack. Noble and Tsien {Jack et al. 1983]). Then. if we assume that only the first 
mode is present, the results don't change qualitatively. 

In order to apply the Haken's slaving principle to Eq. (3.9) with 8, =0 and so to 
obtain Eq. (4.1), previously we must transform this equation in an homogeneous 
one (when 8 =0), i.e. we must transform the additive control parameter (proportio- 
nal to form factor) in a multiplicative one. That is, we must make a change of 
variables that eliminates the inhomogeneous term proportional to form factor. 

When we make this variables change, the new coefficients of linear terms in the 
resulting equations are all negatives, except the coefficient for 4,, that changes of 
sign in the threshold case. Therefore, the slaving principle works in the simple case 
which the only non negligible form factor is F,. 


(a) q, 


Figure 3. The cusp catastrophe in (u, w, q,) space. 
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When the electric current is zero, the equilibrium points of Eqs. (4.2-4.3) are the 
rest state 4,=0 (always stable), the unstable threshold mode amplitude 
4, =[6'-Js]/2c’ and the stable “excited” mode amplitude A, =[b'+Va]/2c'. where 
A=b"4e'(l+ 7X /P). Taking //A, as control parameter, the global bifurcation 
diagram of Eq. (4.2), i.e. Eq. (3.9) for a single mode dynamics. with /(r)=0 . is done 
in Figure 2. We can see there a saddle-node bifurcation, (Guckenheimer & Holmes 
1983; Drazin 1992; Jackson 1989) where d,, corresponds to unstable branch and 
A,, to stable one. 

Now we consider a long constant pulse of intensity {,#0. Then the fiber's 
dynamics changes and new equilibria amplitudes of membrane potential appear, 
that obey a cubic equation with two control parameters, 2, and [,. In fact, 
making the change of variables 


53 
A, =4q, -—— (4.4) 
38 114, 


we obtain the equilibrium equation (where the parameters u and w are functions of 
i, and I.) 


g, +ug,+w=0 (4.5) 


As a consequence of Eq. (4.5) the bifurcation diagram (in codimension two) is 
given now by the folded surface with the cusp catastrophe, (Haken 1978: Jackson 
1989; Britton 1986) defined in (u, w, g,) space, as is shown in Figure 3. 

Using topological and numerical methods, it is possible to extend our analysis 
(see our first paper: Suarez & Sicardi 1996) to cases with N interacting (non slaved) 
modes with non negligible amplitude. The main result is that, depending of the 
influence interval length, there can exist one (stable) or three (two stable and one 
unstable) steady solutions for the Dirichlet’s problem for the N modes. Therefore, 
we can divide the N dimensional mode space in two regions: the decaying set, 
which is the attraction basin of the rest state, and the amplifying set, which is the 
basin of the excited state. The separation surface (see Figure 4) (which we will call 
threshold surface} between the two regions is the stable manifold of the unstable 
steady solution (the threshold steady state). When we apply a current pulse we 
carry the fiber out the rest state. If at the pulse end, the system is in the decaying 
set, the excitation will fail and the fiber returns to the rest point, but if the pulse’s 
end point is in the amplifying set, the system will go to the excited state. Thus. 
each point on the threshold surface corresponds to a well defined threshold 
distribution of membrane polarization along the fiber axis. (We note that while in 
the uniformly polarized case the fiber has a true threshold value of membrane 
potential, in this case it shows a whole family of threshold states of space depen- 
dent membrane polarization. Therefore, for a given fiber, we don't have one single 
threshold, but a family of threshold states. That is, instead of having a well defined 
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30 2s - 
Z=20mm 
Figure 4. The threshold surface for three modes. 


threshold value of membrane potential as in the uniformly polarized case. we have 
a whole family of localized threshold profiles of membrane potential. ) As generally, 
the threshold surface isn’t orthogonal to the dominant mode, the final point can be 
in the decaying set, but with its projection over the dominant mode in the 
amplifying set. So, we can understand in a geometric way, the anodic block from the 
cathode that it is well known experimentally. 

In the following, we return to analysis of full Fitzhugh-Nagumo-Rattay equations 
(3.1-3.2) with a non-null recovery variable. That is, we consider the modal 
equations (3.9-3.10). 


5 - Recovery effects in mode one amplitudes 


We only consider here the evolution equation for the first uncoupled mode. (As 
we note before, for typical values of membrane parameters. we can expect that 
mode one amplitude be significantly greater than the other ones; for a multimode 
analysis see Ref. [Suarez et al. 1997]) Thus. from Eqs. (3.9-3.10), we get the 
following nonhomogeneous equations of the Bonhoeffer-Van der Pol type: (Jack er 
al 1983). 
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ney —— -3,,,4, -aB, +2 w(t), (5.1) . 


— = B(4,-7B,) {5.2) 
at 


When /(1)=0, the horizontal isocline (i. e. the curve such that dB,/dt =0) is given 


by the equation 
¥B, = A, (5.3) 


while the vertical isocline (i. e. the curve such that d4/dt=0) is given by the 
equation 


aB, =-v,A, +69,,,4) - SS a (5.4) 


where 


ni, 
y,=00——"1. (5.5) 
t 


Figure 5. Isoclines for physiological conditions. 


192 


ee ee 


From synergetics to biological excitability 


This nondimensional number depends on the parameters of the electrode-fiber 
system (such as the fiber radius and the electrode -fiber distance). For a suitable 
value of v,, the vertical isocline presents three branches (see Figure 5), with the 
isocline maximum in the upper half plane (8, greater than zero). The middle 
branch connects the left minimum (in 4,,) with this right maximum (in 4,,) of the 
isocline. In this case the system is excitable by a cathodic pulse. 

In physiological conditions, the only rest point (small y) is the origin A, =0. 
8,=0, being in this case the only intersection of the vertical isocline, with the 
horizontal one. However, it is instructive to study first a case with y big enough so 
that there are three intersections of the isoclines (5.3-5.4), that is three rest points. 
two stable focus (one of them is the origin) and one saddle point (see Figure 6). 
Decreasing 7 , the focus different from the origin collides with the saddle point, and 
both disappear (so we reobtain the physiological conditions, see Figure 5). 


40 


Bi 


-20 


o Ala 2 Alb 4 6 


Figure 6. Isoclines for non-physiological conditions. 


Anibal C. Sicardi-Schifino and Roberto E. Suarez Antola 


-2 0 2 4 


Figure 7. The pseudo-separatrix and the vertical isocline. 


In this nonphysiological case (with three rest points), the stable manifold of the 
saddle point is a true threshold separatrix. That is, this curve separates the 
attracting basins of both focus, in the phase plane (4,, 8,). While the origin 
corresponds to the rest state of the fiber, the other focus corresponds to a 
(nonphysiological) true excited state. In this case the behavior of the system is truly 
all or none. When the saddle point and the excited state focus approach first and 
then collide, the separatrix first tends to become parallel to the middle (unstable) 
branch of the vertical isocline (see Figure 7). And then, the separatrix, as such, 
disappears. But, in the orbits behavior (see Figure 8), it remains a memory of the 
previous conditions in the neighborhood of the former separatrix. Even though the 
all or none behavior disappears, something reminiscent of a threshold remains. 
(see the Appendix.) Then, when there is only one intersection point of the isoclines 
(the origin), but y is big enough so that both isoclines are very close in the former 
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Figure 8. The pseudo-separatrix and the orbits behavior. 


intersection region (where the true threshold rest point and the excited state point 
were located in the nonphysiological case), a narrow channel is created and the 
system is attracted to this channel, as it was attracted before to the excited state 
focus. Then the system approaches to the channel from the place of the former 
excited state focus and, after remaining there for a while, it gets out from the place 
of the former threshold side. 

So, this narrow channel works separating the orbits that correspond to full 
action potentials from those which only return to rest. But this channel acts only 
as a fuzzy threshold curve. We call this fuzzy curve. pseudo-separatrix. If we 
consider all the phase space, there is a region where the pseudo-separatix gets 
increasingly fuzzy and finally disappears. However. if we only consider the 
reachable region using electric current pulses from rest. the pseudo-separatrix is 
fairly well defined. So that it can be represented by a curve that locally separates a 
decaying set from an amplifying one. The amplifying set is formed by points that 
originate full action potentials. while the decaying set is formed by points from 
which the system plainly returns to rest. But instead of the decaying and ampli- 
fying sets defined in the case B, =0. these sets only have a local meaning. When the 
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pseudo-separatrix vanishes, the boundary between both sets is fuzzy and the be calculated (Suarez & Sicardi 1996) in terms of a suitable Green function G(r). 
classification of points in decaying or amplifying has not meaning. However, in Thus. the just threshold condition can be formulated as 
physiological conditions, we find again, as in the case B,=0, a whole family of 
threshold states. ' minimum ance —t')dr= 4 =Q, (6.1) 
: te [o.- x) 
1.0 where Q, is the limit threshold charge. The Green function G(t) can be obtained 


using our theory. We show the form of G(r) in two cases (for two different values of 
i 4/1) in Figure 9 and Figure 10. 


same volume, Suarez & Sicardi 1997). This comparison shows a well agreement of 
our theoretical results with experimental observations. Moreover, we have predicted 
new relations between the measurable parameters. Experiments for testing them 


In this way, we can compare our results with experimental ones (see, in this 
0.5 
are in progress. 
1.5 
G(t) 
0.0 
-0.5 | 
' 


Figure 9. The Green function for certain values of parameters. 


6 - Conclusions and biological meaning | 


In clinical electrophysiology the global excitability of a given electrode-tissue-target 
fiber system is usually characterized by the strength-duration curve, i.e. by the 
relation between the duration of an applied rectangular electric current pulse and 
its threshold amplitude. Rheobase, chronaxie (or time constant) and the so called 


© parameter (i.e. the quotient between the threshold amplitude, for pulses whose Time 
duration is equal to system's time constant, and the rheobase) are used as suitable 
Parameters to summarize the experimental results. However. these parameters can Figure £0. The Green function for other values of parameters. 
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In conclusion, we have constructed a modal analysis for the electrode-fiber 
system, in which we take into account the spatial dependence of the external 
electrode field, defining a form factor for each electrode-fiber system. Although here 
we only give the explicit expression for the form factor in the point electrode's case, 
it is possible to calculate the form factor for any type of electrode, using a multipo- 
lar expansion, and in all these cases our theory works. Moreover, in a future work, 
we will extend our theory to other geometries different from a fiber (for example, the 
myocardium). Although here we only consider electric pulses of constant intensity, 
our theory works too for any time dependence of the current intensity /(1). In a 
future work, we will consider a periodic /(¢), in order to study a more interesting 
dynamic behavior, as Hopf bifurcations and the possible emergence of chaotic 
behavior. Also in a future work, we will modify the boundary conditions, in order to 
describe the propagation of the action potential, looking for solutions of soliton-like 
type. 
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Appendix 


We consider in this Appendix the nonphysiological situation of Figure 6. In 
contrast with Figure 5, we have now three rest points for a suitable election of the 


Parameters. When these parameters vary, (for example, if y decreases) we reobtain 
the physiological case with only one rest point. In order to clarify the three rest 
point case, we consider the following mathematical example (where x corresponds 
to 4, and y corresponds to 8,): 


—=-x+3x?-x'-y (Al) 
a (A2) 


We note that the origin is always a rest point. If 8 is positive, the origin is a 
stable focus. When 6 is positive but smaller than 1.25, there are another two rest 
points, which disappear in a reverse saddle-focus bifurcation for B = 125. When B is 
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between 0 and 1.25, one of these points is always a saddle point, while the other 
rest point is first a stable node (for B between 0 and 0.6) and later becomes a stable 
focus (from B=06 on). The stable manifold of the saddle point is the threshold 
separatrix between the two attraction basins of both stable rest points (for 8 < 1.25). 
When § > 1.25, it only exists a single rest point (the origin), but a memory of the 
former situation remains. The reminiscence of the old threshold separatrix is the 
pseudoseparatrix that we consider in the main text. 
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Abstract : For the same nerve or muscle fiber the chronaxy 
measured using a small electrode located inside the fiber, is 
different from the chronaxies measured using the same electrode 
located in the volume conductor outside the fibre. Furthermore, 
measured chronaxies increase with the distance between the 
electrode and the fibre. A mathematical formula that relates the 
time constant of the strength duration curve with the parameters 
of the system composed by the electrode, the volume conductor 
and the fibre is used to explain the aforementioned differences for 
unmyelinated nerve and skeletal muscle fibres. For both 
extacellular and intracellular electrodes, the time constant is an 
increasing function of a characteristic length. This length is 
obtained from a relation between the unstable stationary 
solutions of Nagumo’s nonlinear cable equation and_ the 
polarizing effect of the electrodes on the fibre. In the case of 
extracellular electrodes, the characteristic length is defined as the 
measure of an interval of influence constructed from the 
activating function. In the case of point intracellular 
electrodes,the characteristic length is defined by a construction 
done with the help of the only bounded but not constant 
stationary solution of Nagumo’s equation. 


I. INTRODUCTION 


Let us consider an excitable fibre at its rest state, embedded 
in a volume conductor. An active point electrode is located 
either inside or near the fibre. An indifferent electrode is used 
to close the circuit. In order to produce the emergence of an 
action potential (AP) in the fibre, an electric current pulse can 
be applied to the volume conductor by means of the electrodes. 
For a rectangular current pulse of a given duration, there is a 
threshold amplitude such that a little above it the electric 
stimulation succeeds (an AP is produced), and a little below it 
the stimulation fails. If the duration of the pulse decreases, its 
threshold amplitude (strength) increases, so that the strength- 
duration curve looks like an hyperbola with an horizontal 
asymptote for long pulses and a vertical asymptote for short 
ones [1]. The least threshold current for long pulses is the 
rheobase R. The product of pulse duration (tp) and threshold 
amplitude (Iy) is the threshold charge Qy. This charge 
decreases when tp decreases. There is a limit threshold charge 
(Quo ) for vanishing pulse duration. The time constant is 
defined as 


ts = Quo/R (1) 


The chronaxy tc is defined as the duration of a pulse with a 
threshold amplitude twice the rheobase. tc is closely related 
with ts . Usually both behave in the same way when system’s 
parameters vary. For the same nerve or muscle fibre the 
chronaxy measured using a point active electrode located 
inside the fibre, is different from the chronaxies measured 
using the same electrode located in the volume conductor 
outside the fibre. Furthermore, measured chronaxies increase 
with the distance between the electrode and the fibre [2] [3]. 

In 1975, Rank posed the problem of explaining these 
differences in relation with the electrical stimulation of 
mammalian central nervous system [4]. 

The purpose of the present paper is to give an explanation of 
the aforementioned differences in the chronaxies, for 
unmyelinated nerve fibres and skeletal muscle fibres 
stimulated by point electrodes. 


Il. MATHEMATICAL FORMULAE FOR THE TIME CONSTANTS 


A. Historical background 

Until Hodgkin and Huxley presented their mathematical 
model for the ionic currents in the squid’s giant axon, the 
research on strength-duration curves was mainly experimental. 
After that, digital simulations were possible using 
phenomenological equations for membrane’s ionic currents in 
the framework of non-linear cable theory. 

The first numerical calculations of strength-duration curves 
were done for uniformly polarized membranes and for fibres 
stimulated by internal point electrodes, in the sixties [5]. In 
1976, McNeal simulated a strength-duration curve for a 
myelinated nerve fibre stimulated by an external point 
electrode [6]. In 1986, Rattay showed that the activating 
function that appears in non-linear cable equation is the key 
element to relate the electric current field produced by an 
external electrode with the polarization of fibre’s membrane. 
Using this function, several numerical calculations were done 
for strength-duration curves in extra cellular stimulation [7] 
[8]. 

These numerical works focussed mainly in the rheobase. 
The relation between the time constant (or chronaxy) and the 
parameters of the system composed by the electrodes, the 
tissues and the target fibre, was not studied in any detail. 
Besides, the numerical solution by itself doesn’t yield too 


much insight into the phenomenon of threshold in extracellular 
or intracellular stimulation. 

In 1988, Rubinstein and Spelman published an analytical 
approach to fibre’s polarization produced by external 
electrodes [9]. They used linear cable theory and worked with 
wave numbers in the Fourier domain. They discussed in a half- 
quantitative way some experimental results obtained by 
classical electro physiologists and some more _ recent 
experimental results related with psychophysical thresholds of 
cochlear implants. They introduced an ad-hoc assumption 
about threshold, because a linear approach can not cope with 
thresholds and its relations with systems parameters. 

In 1994 appeared a fully non-linear analytical approach to 
threshold phenomena in biological tissues stimulated by 
external electrodes [10] [11]. In the fibre’s case an unbounded 
spatial domain problem is transformed into a bounded domain 
problem, introducing the concept of interval of influence of the 
electrode on the fibre. This allows the application of methods 
of non-linear modal analysis closely related with Synergetic 
[12]. 


B. Theoretical results for time constants obtained from non- 
linear modal analysis of the stimulation process by 
extracellular electrodes 

In order to simplify as much as possible the mathematics of 
threshold phenomena, let us consider a fibre in which 
membrane’s recovery variables are frozen at their rest values, 
but membrane’s activation variables are always fully relaxed 
to equilibrium with the local value of membrane’s voltage 

(Vm). Then the density of ionic current through the membrane 

is only a function J(Vy) of membrane voltage. A cubic 

polynomial with three real roots (Vp , Vy and Vg) may be a 

good approximation. Vp is the rest state (stable), Vy is the 

threshold for uniform membrane depolarization (unstable), and 

Vz is an excited state (stable, due to the absence of recovery 

effects). 

In this case the non-linear cable equation reduces to the 

following Nagumo-Rattay equation [11]: 


0°V 


ox? 


oV 
T™ ° a +Ry -J(VM) =A { 
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Here ty = RyCqm and A°y = (Ry 0, ap)/2 . Cy is the 
capacitance per unit area of membrane, Ry is a resistance per 
unit area of membrane defined by 


1 _ dt(Ve) 3) 
Ry dV 


Then RyJ(Vq)=(Vy-Vp)-b(V-Vp) +¢(Vy-Vp)* , where b and 
c are suitable positive parameters. ©, is the conductivity of the 
intracellular space and ag is the radius of the fibre. F(x) i(t) is 
Rattay’s activating function, being i(t) the applied electric 
current pulse and F(x) a function of the coordinate x along the 
fibre. Fig.1 shows the construction of an interval of influence 
of the electrode on the fibre. In order to study threshold 


dynamics, F(x) may be neglected outside this interval. [10] 
[11] [13]. At the end points of that interval, membrane voltage 
can be considered as fixed at its rest value. So, Dirichlet’s 
boundary conditions can be applied to the field Vy(t,x), during 
membrane polarization, until a threshold voltage pattern is 
attained. 


Fig.1. Construction of the interval of influence using the activating function. 


Reference [11] shows that for cathodic make stimulation by an 
external electrode, the following formula for the time constant 
is obtained applying non-linear modal analysis to Eq. (2) : 


on ae (4) 


Here ¢ is the length of the interval of influence. When the 


volume conductor is ohmic, homogeneous and isotropic ¢ =6z 
, being z the distance between the point electrode and the fiber. 
Substituting ?=6z in Equation (4) it follows that : 


6 2 

ty ( | 

= 5 (5) 
M 14(6%/) 


ts is an increasing function of the electrode-fiber distance. 
When z grows without bounds, ts tends to Ty. When z 
approaches zero, ts tends to zero (Fig. 2). 


x 


Fig.2. Relation between Y = ts/ Ty and X=6z/TAy 


C. Theoretical results for time constants measured with 
intracellular electrodes 

Equation (4) was obtained for extracellular electrodes. 
However, if the parameter @ is suitably redefined, this formula 
can be applied to the strength-duration curves produced by 
intracellular point electrodes. A clue to find the new meaning 
of & is obtained from the stationary patterns of polarization 
corresponding to a fibre described by Nagumo‘s non-linear 
cable equation (that is, by Equation (2) without the activating 
function). In the case of extra cellular stimulation, for each 
interval of influence there is a stationary threshold pattern of 
membrane polarization. To this threshold pattern corresponds a 
threshold amplitude of the dominant mode that results from 
non-linear modal analysis. In order to produce an AP, the 
polarization due to the external cathode must leave the 
dominant mode amplitude above its threshold. The 
abovementioned threshold patterns are unstable stationary 
solutions of Nagumo’s equation. In all cases there is an 
interval centred in the point of the fibre that is nearest to the 
electrode, such that there the local values of membrane voltage 
are above the threshold for uniform membrane depolarization 
(Vy). As a consequence, a necessary condition that a voltage 
pattern must fulfil to be a threshold one is that a liminal length 
of fibre membrane has Vy above Vy. Thus, the concept of 
liminal length for nerve or muscle excitation by external 
electrodes, introduced in 1937 by Rushton, is obtained as a by- 
product of this analysis. 

When the stimulation is produced by an internal anode, 
there is a single stationary threshold pattern of depolarization. 
The pattern is symmetric, centred in the point of the fibre in 
which the electrode is located and is fairly well localized. As 
shown in Fig. 3, using this threshold pattern, an interval of 
length ?, can be constructed. At the end points of this interval, 


it is possible to establish that Vy=Vp .This assumption can be 
done during the depolarization of the membrane produced by a 
current pulse below or just above threshold (before an AP is 
produced). A non-linear modal analysis can be done. 
Threshold amplitude for a dominant mode is obtained, which 
is closely related with the projection of the threshold pattern 
onto the corresponding eigenfunction. Working as in the case 
of extra cellular stimulation, Equation (4) for ts is obtained 


with @, in place of £ | In order to obtain an estimation of ¢ 


an expression for the stationary threshold pattern for 
intracellular stimulation must be considered. A_ useful 
approximation to this pattern was obtained by D. Noble in 
1972 [5]. 


-(xy +Am) 


Fig.3. Construction of an interval that gives the characteristic length for a 
point intracellular electrode. 


Using this approximation, and following the construction 
procedure suggested in Fig. 3, an analytical estimation is 
obtained: 


e; = 2(xy +A) (6) 


Here 2xy is Rushton’s liminal length as defined by D. 
Noble. 
Substituting Equation (6) in Equation (4) in the place of ¢ , a 
formula for the time constant for intracellular point electrodes 
is obtained: 


T 
t,=— (7) 
2 
Tt 1 


[eX J} 


ts is an increasing function of xy/Ay . But, as D. Noble has 
shown [5] 


= (8) 


As a consequence of Eq. (7) and Eq. (8), ts is a function of the 
slope of the density of membrane’s ionic current when 
Vu=Vy. When this slope increases (decreases) in absolute 
value, ts decreases (increases). 


Il. CONCLUSIONS 


A. Comparison between the time constants for a fiber 
stimulated by extracellular and intracellular electrodes 
Eq. (4) applies equally to excitation by external and by 
internal electrodes. It follows that when the length of the 
interval of influence of a point external electrode ¢= 6z is 


equal to @ ; for the stationary threshold pattern corresponding 


to a point internal electrode, both time constants will have the 
same value. So, when the distance between the point external 


electrode and the fiber verifies zy = (1/6) @ ; » both chronaxies 


will be equal or approximately equal. We obtain from 
Equation (6) the following expression for this distance Zo : 


1 


When z is greater (less) than Zo, the time constant for a fibre 
stimulated by a point external electrode will be greater (less) 
than the time constant of the same fibre stimulated by a point 
internal electrode. This result may explain the experimental 
results about the time constant (or chronaxy) mentioned in the 
introduction of this paper. Moreover, we have a quantitative 
prediction that could be tested by doing suitable experiments 
or digital simulations. 

B. Time constants for intracellular stimulation 

For the giant axon of the squid: xy/Aq ~ 0.8. Then from 
Equation (6) it follows that ts/t = 0.57. In this case ty = 
3mseg, so that the estimated value of ts is 1.7mseg. . The 
measured value is 2.0mseg. . The estimation is fairly good, 
taking into account that the model for membrane’s ionic 
current used is highly simplified. 

For a typical skeletal muscle fibre, xy/Ay, = 0.4, so that ts/Tyy 
= 0.44. For a Purkinje fibre of the heart’s conduction system, 
Xy/Am = 0.15, so that ts/ty = 0.35. In all these cases the 
theoretical relation between ts/ty and xy/Ay gives values fairly 
near the experimentally determined ones. However, the error 
between the estimation and the measurement seems to increase 
when Xy/Ay decreases [10]. 


C. Possible improvements 
In the theoretical expression for ts the estimation of ¢; may 


be improved substituting Ay, by a somewhat greater space 


constant nie This would allow us to take into account the 


effect of the nonlinearities in the stationary threshold pattern 
for internal electrodes. This is not taken into account in the 
analytical approximation proposed by D. Noble. The modified 
expression for ts would give better estimations for the time 
constants. 

To apply the modal approach to the experiments done by 
classical electro physiologists, it is necessary to modify 
Nagumo-Rattay’s equation because in this case the fibre is 
wrapped by a sheath of connective tissue with a conductivity at 


least an order of magnitude smaller than the conductivity of 
the volume conductor (Ringer solution) [10] [14]. 

It is possible to apply the methods of non-linear modal 
analysis to a cable model with an ionic current density with 
free activation and recovery variables, obtaining thus an 
improved theoretical expression for time constants, both for 
internal and external electrodes [10] [13]. 


Then, a similar comparison as was done in this paper could 
be made. 
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Abstract : In the framework of cardiac pacing, measured 
chronaxies are different for current and for voltage pulses. 
They increase with the distance between the pacing cathode 
and the myocardium, and with the size of the electrode. 
Mathematical formulae that relate the time constants (and the 
chronaxies) of the strength duration curves with the 
parameters of the system composed by the electrode, the 
volume conductor and the myocardium are used to explain the 
aforementioned experimental results and to discuss, from an 
engineering standpoint, some aspects related with thresholds 
for cardiac pacing, both for unipolar and bipolar pacing. These 
formulae are derived from a mathematical model of the electric 
stimulation of a syncytium framed in the bidomain theory. The 
time constants are given as increasing functions of a non 
dimensional number, equal to the quotient between a a certain 
characteristic length and myocardium’s space constant. In the 
case of pacing electrodes, the characteristic length is a function 
of the size and shape of a region of influence of the cathode on 
the excitable tissue located adjacent to the electrode. 


I. INTRODUCTION 


The myocardium is composed by a mass of excitable cells 
tightly coupled, so that an electrical syncytium is formed [1]. 
Let us consider an active cathode located in a volume 
conductor, near the mass of excitable cells that make the 
tissue (an active electrode at the tip of a lead for cardiac 
pacing, for example). An anode is used to close the circuit (a 
ring electrode in the lead mentioned before, perhaps at 2cm 
or less from the tip, or a truly indifferent electrode located in 
the pulse generator, perhaps at 5 to 10cm from the tip) [2]. 
For a rectangular cathodic current pulse of a given duration, 
there is a threshold amplitude such that a little above it the 
electric stimulation succeeds (an AP is produced), and a little 
below it the stimulation fails. If the duration of the pulse 
decreases, its threshold amplitude (strength) increases, so 
that the strength-duration curve looks like a hyperbola with a 
horizontal asymptote for long pulses and a_ vertical 
asymptote for short ones[3]. The results of many laboratory 
experiments and clinical practice show that the parameters of 
the strength-duration curve depend not only of the 
excitability properties of the membranes of the target 
element (myocardium), but depend also of certain geometric 
and bioelectric parameters of the system formed by the 
electrodes, the tissues and the target element[4],[5]. The 
least threshold current for long pulses is the rheobase /,p, 


and the product of pulse duration (f,) and threshold 


amplitude (/;) is the threshold chargeQ,;. This charge 


decreases when tp decreases. There is a limit threshold 
charge (Q;.) ) for vanishing pulse duration. Then the time 


constant for electric current pulses is defined as 


te, ote 
SI i 

The chronaxy tc is defined as the duration of a pulse with 
a threshold amplitude twice the rheobase. It results that tc is 
closely related with ts. Usually both behave in the same way 
when system’s parameters vary. For the same myocardium, 
the chronaxy measured using a point active electrode, acting 
as anode and located inside a cell, is different from the 
chronaxies measured using the same electrode acting as 
cathode and located in the volume conductor outside the 
tissue. Furthermore, for external electrodes, measured 
chronaxies increase with the distance between the electrode 
and the tissue and with electrode’s size. When a pacing 
electrode is located near enough to myocardium, measured 
chronaxies depend of electrode’s shape [6]. 

In cardiac pacing for medical purposes, a rectangular 
voltage pulse is applied to the electrodes, instead of the 
rectangular current pulse mentioned before [4], [7]. Due 
mainly to the polarization of the electrode-tissue interface, 
the electric load seen by the pulse generator is not an ohmic 
one. The amplitude of the current pulse now decreases from 
the beginning to the end of the pulse while the amplitude of 
the voltage pulse remains constant. A strength duration 
curve may be measured for threshold voltage pulses, and a 
voltage rheobaseUp,, a threshold value of the product 


between pulse duration and_ threshold voltage 
amplitude U (-,). and a voltage chronaxy (or a time 


(1) 


constant ft, ;; ) may be determined from this measured curve. 


The purpose of this paper is twofold: (a) to extend a 
previous work reported in reference [8] on the time constants 
for cathodic make excitation of nerve and skeletal muscle 
fibres, to embrace an electric syncytium like myocardium, 
and (b) to discuss the relation between thresholds for 
rectangular current and voltage pulses from an engineering 
standpoint. 

The approach to the strength-duration curves by people 
working in cardiac pacing has been mainly an empirical one. 
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The intended theoretical derivations of well known empirical 
formulae like Lapicque’s (1907) or Weiss’ (1901) for 
myocardial stimulation by pacing electrodes were sometimes 
very ingenious, as the one presented in reference[7], but not 
in accord with standard electrophysiological knowledge. In 
1977 Lindemans obtained experimental data that suggest 
that an initial depolarization, above uniform membrane 
threshold, of a critical mass of cardiac cells is needed to 
produce a propagated AP [9]. This pointed towards the 
existence of liminal volumes for the excitation of a three 
dimensional syncytium, analogous to the liminal lengths 
introduced by Rushton (1937) in the case of nerve or muscle 
fibres [8]. However, it lacked an appropriate theory to give 
an explanation of these findings. Under the name of 
bidomain theory of electro-cardiology, in 1978 began the 
construction of a sound, albeit phenomenological theoretical 
framework [10]. Most work was done in applying bidomain 
theory to digital simulation of AP propagation and of cardiac 
far-field stimulation (due to its importance in cardioversion 
and defibrillation), and comparatively much less effort was 
done on near-field stimulation, which is the situation of 
interest for cardiac pacing [11], [12]. In 1994 appeared a 
fully analytical non-linear approach to threshold phenomena 
in biological tissues stimulated by external electrodes [8], 
[13], [14], [15]. This approach, framed in the bidomain 
model, allows an explanation of the liminal volumes 
discovered by Lindemans, as well as a derivation of their 
geometric dimensions. The polarizing effect of an external 
electrode on the excitable electrical syncytium is described 
by a bounded region of influence. This allows the 
application of methods of non-linear modal analysis closely 
related with the methods of Synergetics [16]. 


Il. MATHEMATICAL FORMULAE FOR THE TIME 
CONSTANTS 


A Time constants for current pulses. 

In order to simplify as much as possible the 
mathematics of threshold phenomena, let us consider first a 
homogeneous and isotropic electric syncytium. Let us 
further assume that membrane’s recovery variables are 
frozen at their rest values, but membrane’s activation 
variables are always fully relaxed to equilibrium with the 
local value of membrane’s voltage. Then the density of ionic 
current through the membrane is only a function of 
membrane voltage. A cubic polynomial, with three real 
roots, may be a good approximation for the cardiac muscle 
[17]. In this case the three dimensional non-linear Cable 
Equation of the continuum of excitable membranes reduces 
to the Nagumo Equation. Assuming that in the bidomain’s 
boundary the electric current goes to or comes from the 
surrounding tissues only through the interstitial domain, in 
the surface of the cardiac muscle the current density in the 
intracellular continuum must vanish. The region of influence 
of the applied current field over the membranes of the 
myocardium cells may be represented by half an ellipsoid of 
revolution [13], [15] (Fig1). 


Fig.1. Electrode, myocardium and region of influence 


The radius R, of the circle located on the boundary of the 
excitable tissue, nearest to the cathode, is given by 

Ro = (1 + e) (2) 
7% is a cathode’s radius, € is the thickness of the layer of 
non-excitable tissues interposed between the electrode and 
the myocardium, and y is a non-dimensional parameter 
(related with the spatial distribution of the electric current 
density) that will be explained below. The minor axis of the 
half-ellipsoid, cy, that gives a measure of the depth of 
penetration of the polarizing effect in the syncytium, is a 
function of Ry and of the syncytium’s space constant A, 
[13], [15]: 

Co B(Ro Au ) 


Am a,” +(Ry/Ay 


(3) 


(Here # is near 1). ~ 5.783. It can be seen thatcyis a 
monotonically increasing function that tends to £.A,, when 
the radius of the half ellipsoid of influence goes to ©. A 
necessary condition for AP generation is that the region of 
influence of the electrode includes a liminal volume of 
syncytium. SoRy >Ro,., being Roa minimum value for 
excitation [13], [15]. Membrane voltage is supposed to be at 
its rest value on the curved surface of the half-ellipsoid, and 
the normal derivative of membrane voltage on the flat face is 
proportional to the applied cathodic current J (r). Using this 
region of influence it is possible to carry on with the 
methods of non-linear modal analysis, and after doing that 
and making several simplifications, the following results are 
obtained [13], [15]. 

A current pulse of duration f,, is just threshold, 1/7 (t) , if 


j ol p -#)I r(thdt = Oro (4) 


t 


Qry is the limit threshold charge, G(t)=e ‘ (5) 
ey ae) 


: 6 
oe a+blRy2/Ay?) ) 


met 
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Ty is the linear membrane time constant with the activation 
variables relaxed to equilibrium with membrane voltage, Rp 
and A, were defined before, and the non-dimensional 
parameters a and b may be takena = 38.41 and b ~ 6. If the 
amplitude of the current pulse is constant/;,9, from 


Equations (4) and (5) we obtain Lapicque’s formula for 
cathodic make excitation of an electric syncytium, like 
myocardium: 

I 
Iroltp)=—4- o 


'p 


1- e's 
The rheobase is Ip = so that ¢; , given by formula 
Sl 


(6) as a well defined function of the parameters of the 
electrode, of the excitable tissue, and of the geometry of the 
electric current field, is the time constant; , for cathodic 


make excitation. Figure 2 shows the relation between the 


t 
dimensionless variables A and Ry, : 
Tu Au 


Fig.2.Relation between Y =f, , Jee and X =R,/Ay 


The time constant increases with the radius of the 
electrode, with the thickness of the layer of non-excitable 
tissues and with the parameter vy . For unipolar stimulation by 


a convex cathode,y~1.33, but it increases in bipolar 


stimulation when the distance between the ring anode and 
the cathode decreases [13], [15]. All this is in accord with 
experimental results [4], [7], [18], [19], [20]. 


B Time constants for voltage pulses. 

Modern permanent pulse generators use constant voltage 
pulses to pace the heart. While the voltage remains at the 
programmed value, the electric current suffers variations that 
reflect the changes in the load impedance seen by the pulse 
generator. This impedance is the sum of an ohmic 
resistance R,,and the impedance produced by _ the 


polarization of the electrochemical double layer located at 
the interface between the electrodes and the tissues [4], [7], 
[18]. The voltage drop between the anode and the cathode 
U (r) can be written, in a linear approximation, as the ohmic 
voltage drop in the cables and the tissues plus the voltage 
drop in the double layer: 


u(t)=R, 1(t)+ | K(t—u)t(w)du (8a) 


K(t)is positive and [Kudu =R,, R,is the D.C. 
0 


resistance of the interface. As U (t) is given by the voltage 


pulse generator, it is convenient to invert (8a), thus 
obtaining: 


i= (v0) [ul “ule (8b) 
R, 0 
a6 R 
L(t) is also positive and [ L(u)du =— . Fora 
0 


R ae Ry 
voltage pulse of constant amplitude Uy, Equation (8b) may 
be rewritten thus: 


Uo 
I(t)=——— 1 
(") one (1+ (7) (9) 
; o ; ’ Ra 
p(t) is positive, monotonically decreasing from ~(0)= “a 


Pp 


towards zero when f— oo. It can be given as a function 


t 
of [ L(u)du . From (9) it follows that J (r) decreases from 
0 


the beginning to the end of the voltage pulse. A voltage 
pulse of constant amplitude will be a threshold pulse U; 9 if 
the corresponding current pulse J, (r) verifies Equation (4). 
After having substituted the threshold current pulse given by 
Equation (9) forU; 9, in Equation (4), we 


obtain: 


Urolt,)=(Rp + Ry }Oro/({ (art {ale, —2ole)an (10) 
0 0 


t 


IfG(t)=e it is possible to show __ that: 


ty - 
lim, = [al, ~rlo(t)at =0. But i G(t)at =t, 7, So taking 
0 0 


into account that Ip ak , the rheobase of voltage 
Si 


is: Up =lim, _,.. Urolt, )=(R, +a Ve (11) 
Voltage rheobase and current rheobase are connected 

by the total resistance (D.C. impedance) seen by the pulse 

generator. 

The strength-duration formula for rectangular voltage 

pulses, given by Equation (10), may be rewritten thus: 


Dus (.,)= ts,UpR 


fates fo, =A 


The strength-duration formula for rectangular current 
pulses obtained from Equation (4) may be rewritten thus: 


(12) 
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troy) = (3222) (13) 


ty 


J G(t).dt 


For short pulses they both behave proportional to iS , but for 


longer pulses the threshold amplitude of rectangular voltage 

pulses approaches to its rheobase faster than the threshold 

amplitude of rectangular current pulses of the same 

durations. All this is in accord with experimental results [7], 

[18], [20]. 

From Equation (12), taking into account that, by definition, 
lim, _,0 (t, Ur oltp 


tsy = , it follows that: 
UR 


ts 7 R 


P 
a 1+9(0) R,+R, st ”) 

Thus the time constant for excitation with voltage 
pulses is always less than the time constant for excitation 
by current pulses, and the difference between them 
increases with the increase in the D.C. impedance of the 
electrochemical double layer relative to the ohmic 
impedance of the lead cables and the tissues. All this is in 
accord with known experimental results [4], [7],[18], [19], 
[20] (Figure 3) 


nme sil 


7 


mm lo 


Fig.3. Relation between ts and ro (modified from [7]). 


Ill. CONCLUSIONS 

The theoretical results for the time constants, both for 
current pulses and voltage pulses agree fairly well with 
known experimental results from  electrophysiological 
experiments and from clinical practice. The quantitative 
predictions implicit in the formulae given in this paper 
should be tested against laboratory experiments or digital 
simulations of the emergence of an AP in electrical syncytia 
like cardiac tissue. Also, they could help in the design both 
of complex digital simulations of threshold phenomena and 
in the assessment of electrodes for chronic cardiac pacing. 

The mathematical model was constructed for a 
homogeneous and isotropic syncytium. As the cardiac 
muscle is heterogeneous, anisotropic and with unequal 
anisotropy (the anisotropy ratios of the intracellular and the 
interstitial continua differ) the analytical approach developed 
here should be extended to cope with this anisotropy. 


The time scales of the pulses used in cardiac pacing justify 
the assumptions made about the relaxation of the activation 
variables to equilibrium with membrane voltage and the 
freezing of the recovery variables. However, it is possible to 
take into account activation and recovery effects in the 
framework of a nonlinear modal analysis of threshold 
phenomena [13], [15], [16]. 
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